set.seed(12618)
Ejemplo ilustrativo de la verosimilitud de una Bernoulli vs nuestra Bernoulli truncada
Graficar el sesgo del EMV en función del tamaño muestral n suponiendo que θ = 1/4.
tita <- 1/4
p_real <- 1/6 + 2/3 * tita
Nrep <- 10000
n_vals <- c(1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 100, 200, 300, 400, 500, 1000, 1500, 2000)
sesgos <- numeric(length(n_vals))
esperanza_p <- function(n, p_real, k_esimo_momento = 1) {
si_promedio_menor_o_igual_a_un_sexto <- function(n, p_real, k_esimo_momento = 1) {
k <- floor(n / 6)
return(((1/6)**k_esimo_momento) * pbinom(k, size = n, prob = p_real))
}
si_promedio_entre_un_sexto_y_cinco_sextos <- function(n, p_real, k_esimo_momento = 1) {
lower <- floor(n / 6) + 1
upper <- ceiling(5 * n / 6) - 1
if (lower > upper) return(0)
k_vals <- lower:upper
probs <- dbinom(k_vals, size = n, prob = p_real)
pesos <- (k_vals / n) ** k_esimo_momento
return(sum(probs * pesos))
}
si_promedio_mayor_o_igual_a_cinco_sextos <- function(n, p_real, k_esimo_momento = 1) {
lower <- ceiling(5 * n / 6)
return(((5/6)**k_esimo_momento) * (1 - pbinom(lower - 1, size = n, prob = p_real)))
}
return(
si_promedio_menor_o_igual_a_un_sexto(n, p_real, k_esimo_momento) +
si_promedio_entre_un_sexto_y_cinco_sextos(n, p_real, k_esimo_momento) +
si_promedio_mayor_o_igual_a_cinco_sextos(n, p_real, k_esimo_momento)
)
}
esperanza_tita <- function(n, p_real) {
return ((3/2) * (esperanza_p(n, p_real, 1) - 1/6))
}
for (i in seq_along(n_vals)) {
n <- n_vals[i]
estimadores <- numeric(Nrep)
sesgos[i] <- esperanza_tita(n, p_real) - tita
}
El gráfico muestra el del estimador de máxima verosimilitud para distintos tamaños muestrales \(n\), con \(\theta = \frac{1}{4}\). Es más pronunciado cuando \(n\) es pequeño. A medida que \(n\) aumenta, el estimador tiende a concentrarse en torno a su valor esperado y el sesgo disminuye, tendiendo a cero.
Del gráfico podríamos deducir que el estimador resulta ser asintóticamente insesgado.
Graficar el ECM del EMV en función del tamaño muestral n suponiendo que θ = 1/4. Superponer un gráfico del ECM del EMV calculado utilizando el método clásico suponiendo que la proporción de estudiantes que se copió alguna vez en un examen y miente cuando se le hace la pregunta directa es π = 1/3.
proba_pi <- 1/3
p_y_clas <- (1 - proba_pi) * tita
Nrep <- 10000
n_vals <- seq(2, 200, by = 2)
ecm_aleatorizado <- numeric(length(n_vals))
ecm_clas <- numeric(length(n_vals))
ecm_aleat <- function(n, p_real){
sesgo_tita <- function(n, p_real){
return (esperanza_tita(n, p_real) - tita)
}
varianza_tita <- function(n, p_real){
return (((3/2)**2) * (esperanza_p(n, p_real, k_esimo_momento = 2) - esperanza_p(n, p_real)**2))
}
return (varianza_tita(n, p_real) + sesgo_tita(n, p_real)**2)
}
for (i in seq_along(n_vals)) {
n <- n_vals[i]
ecm_aleatorizado[i] <- ecm_aleat(n, p_real)
muestras_clas <- rbinom(Nrep, n, p_y_clas)
promedios <- muestras_clas / n
ecm_clas[i] <- mean((promedios - tita)^2)
}
El gráfico compara el del estimador de máxima verosimilitud (EMV) obtenido mediante el método de (curva azul), con el ECM del estimador (curva roja), suponiendo que la proporción verdadera de individuos con la característica de interés es \(\theta = \frac{1}{4}\) y que el porcentaje de estudiantes que mienten en el método clásico es \(\pi = \frac{1}{3}\).
Método numérico
p <- seq(1/6, 5/6, length.out = 20000)
l_diff <- -2 * 20 * log(p) - 2 * 80 * log(1 - p) +
2 * 20 * log(20 / 100) + 2 * 80 * log(80 / 100)
tita <- 3/2 * (p - 1/6)
corte <- qchisq(1-0.05, df = 1)
Simular ambos experimentos en R suponiendo θ = 1/4 y p = 1/3, es decir, generar las respuestas de los encuestados por el método clásico y por el método de respuesta aleatorizada, para n = 10, n = 100 y n = 1000 y calcular el ECM empírico. Presentar los resultados en una tabla y describirlos, verificando que son compatibles con los resultados teóricos hallados en los ítems anteriores.
tita <- 1/4
proba_pi <- 1/3
p_y_aleat <- 1/6 + 2/3 * tita
p_y_clas <- (1 - proba_pi) * tita
n_vals <- c(10,100,1000,2000,3000)
Nrep <- 10000
ecm_aleat <- numeric(length(n_vals))
ecm_clas <- numeric(length(n_vals))
for (i in seq_along(n_vals)) {
n <- n_vals[i]
estim_aleat <- numeric(Nrep)
estim_clas <- numeric(Nrep)
for (iter in 1:Nrep) {
muestra_aleat <- rbinom(n, 1, p_y_aleat)
media_aleat <- mean(muestra_aleat)
estim_aleat[iter] <- min(1, max(0, 1.5 * (media_aleat - 1/6)))
muestra_clas <- rbinom(n, 1, p_y_clas)
estim_clas[iter] <- mean(muestra_clas)
}
ecm_aleat[i] <- mean((estim_aleat - tita)^2)
ecm_clas[i] <- mean((estim_clas - tita)^2)
}
tabla_resultados <- data.frame(
n = n_vals,
ECM_Clasico = ecm_clas,
ECM_Aleat = ecm_aleat
)
print(tabla_resultados)
## n ECM_Clasico ECM_Aleat
## 1 10 0.020782000 0.0415575000
## 2 100 0.008346730 0.0050118025
## 3 1000 0.007087487 0.0005029204
## 4 2000 0.006994958 0.0002563438
## 5 3000 0.007003293 0.0001653930
B <- 100000
y <- c(rep(1, 20), rep(0, 80))
alpha <- 0.05
z_alpha_2 <- qnorm(1 - alpha/2)
theta_hat <- function(sample) {
p_hat <- mean(sample)
theta <- min(1, max(0, 3/2 * (p_hat - 1/6)))
return(theta)
}
theta_bootstrap <- replicate(B, {
sample_y <- sample(y, size = 100, replace = TRUE)
theta_hat(sample_y)
})
IC_bootstrap <- quantile(theta_bootstrap, probs = c(alpha/2, 1 - alpha/2))
sd_boot <- sd(theta_bootstrap)
cota_inf <- theta_hat(y) - z_alpha_2 * sd_boot
cota_sup <- theta_hat(y) + z_alpha_2 * sd_boot
IC_bootstrap_2 <- c(max(0,cota_inf), min(1,cota_sup))
## 2.5% 97.5%
## 0.00 0.17
## [1] 0.00000 0.14908
calcular_nivel_empirico <- function(theta0 = 1/4,
alpha = 0.05,
B = 10000,
n_vals = c(200, 300, 400, 500, 1000, 1500, 2000)) {
crit_val <- qchisq(1 - alpha, df = 1)
p0 <- 1/6 + (2/3) * theta0
computar_T <- function(y) {
n <- length(y)
S <- sum(y)
ybar <- S / n
p_hat <- min(5/6, max(1/6, ybar))
T <- - 2 * (S * log(1/3) + (n - S) * log(2/3) - S * log(p_hat) - (n - S) * (log(1 - p_hat)))
return(T)
}
niveles_empiricos <- numeric(length(n_vals))
for (i in seq_along(n_vals)) {
n <- n_vals[i]
# B muestras de tamaño n bajo H0
muestras <- matrix(rbinom(n * B, size = 1, prob = p0), nrow = B, ncol = n)
T_vals <- apply(muestras, 1, computar_T)
niveles_empiricos[i] <- mean(T_vals > crit_val)
}
return(data.frame(n = n_vals, nivel_empirico = niveles_empiricos))
}
# n chicos
n_pequenos <- 1:150
resultado_pequenos <- calcular_nivel_empirico(n_vals = n_pequenos)
# n grandes
n_grandes <- c(200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000, 2500, 3000)
resultado_grandes <- calcular_nivel_empirico(n_vals = n_grandes)
n_comb <- c(n_pequenos, n_grandes)
nivel_comb <- calcular_nivel_empirico(n_vals = n_comb)
En estos gráfico se observa que:
Este resultado muestra que, si bien el test tiene un comportamiento deficiente para muestras pequeñas, su nivel de significación se ajusta al valor teórico a medida que aumenta el tamaño muestral. Por lo tanto, es recomendable utilizar este test con \(n\) suficientemente grande.
n <- 100
B <- 10000
alpha <- 0.05
theta0 <- 1/4
p0 <- 1/6 + (2/3) * (theta0)
c_crit <- qchisq(1 - alpha, df = 1)
p_vals <- sort(c(seq(1/6, 5/6, by = 0.01), 1/3))
potencia <- numeric(length(p_vals))
calc_aux_logver <- function(p, ybar) {
ybar * log(p) + (1 - ybar) * log(1 - p)
}
for (i in seq_along(p_vals)) {
p <- p_vals[i]
stats <- replicate(B, {
y <- rbinom(n, 1, p)
ybar <- mean(y)
p_hat <- min(5/6, max(1/6, ybar))
stat <- - 2 * n * (calc_aux_logver(p0, ybar) - calc_aux_logver(p_hat, ybar))
stat
})
potencia[i] <- mean(stats > c_crit)
}